Water structure-forming capabilities are temperature shifted for different models 
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A large number ol water models exists for molecular simulations. They differ in the ability to 
reproduce specific features of real water instead of others, like the correct temperature for the density 
maximum or the diffusion coefficient. Past analysis mostly concentrated on ensemble quantities, 
while few data was reported on the different microscopic behavior. Here, we compare seven widely 
used classical water models (SPC, SPC/E, TIP3P, TIP4P, TIP4P-Ew, TIP4P/2005 and TIP5P) in 
terms of their local structure-forming capabilities through hydrogen bonds for temperatures ranging 
from 210 K to 350 K by the introduction of a set of order parameters taking into account the 
configuration of the second solvation shell. We found that all models share the same structural 
pattern up to a temperature shift. When this shift is applied, all models overlap onto a master 
curve. Interestingly, increased stabilization of fully coordinated structures extending to at least 
two solvation shells is found for models that are able to reproduce the correct position of the 
density maximum. Our results provide a self-consistent atomic-level structural comparison protocol, 
which can be of help in elucidating the influence of different water models on protein structure and 
dynamics. 



I. INTRODUCTION 

At the fundamental level, water directly influences sev- 
eral biologically relevant processes including protein fold- 
ing pp, protein-protein association PHI] and amyloid ag- 
gregation [5]. Surprisingly, relatively simple models with 
fixed charges and geometry are able to reproduce the 
phase diagram as well as many of the anomalies of wa- 
ter with good accuracy [H [7] . For example, all popular 
classical water models present a density maximum [5J [3] . 
However, only those that explicitly included this informa- 
tion in the fitting of the potential are able to correctly re- 
produce the experimental value located at around 277 K 
at ambient pressure [10] . 

Due to their improved speed, biomolecular simulations 
in explicit water were traditionally run with TIP3P[TTj or 
SPC [12]. Nowadays, more elaborated models can be eas- 
ily used and their impact on the calculation assessed [13] • 
Optimized four site models reproducing the experimen- 
tal temperature of maximum density seem to improve 
the quality of biomolecular simulations. Best and collab- 
orators showed that predicted helical propensities are in 
better agreement with experiments when a TIP4P/2005 
water model is chosen in place of the traditional TIP3P 
|14j . Others reported that TIP4P-Ew provides better 
free-energy estimations compared to conventional water 
models [15] . In both studies, the improved behavior was 
not connected to a clear microscopic property of the wa- 
ter model. To this aim, one limitation is the lack of a 
common framework to compare the structural behavior 
of liquid water at the atomic level. 

Here, seven popular classical water models, namely 
SPC [ig, SPC/E PI], TIP3P pi], TIP4P pi], TIP4P- 
Ew [IFJ, TIP4P/2005 [5] and TIP5P PI] are investi- 
gated in terms of their local structure forming capabili- 
ties. That is, their ability to form structured or partially 
structured environments of the size of up to two solvation 
shells through hydrogen bonds. This approach represents 
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FIG. 1. Schematic representation of the fully coordinated 
water configuration (P4 population, see text). Dashed lines 
represent hydrogen bonds 

a simplified version of the recent complex networks analy- 
sis introduced to study water structural inhomogeneities 

eoi mi- 
ll. METHODS 

A. Simulations 

Molecular dynamics simulations were run with the pro- 
gram GROMACS [22] with an integration time-step of 
2 fs. The water box consisted of 1024 water molecules 
in the NPT ensemble with pressure of 1 atm and tem- 
peratures ranging from 210 K to 350 K with steps of 10 
K. The Berendsen barostat [23J , velocity rescale thermo- 
stat [21] and PME [22] were used for pressure coupling, 
temperature coupling and long-range electrostatics, re- 
spectively. 

The present analysis was done over 25'000 snapshots 
obtained from a 100 ps long run after a 10 ns equilibration 
in the same conditions. The short length of the analyzed 
trajectory is justified because structural properties are 
calculated for each water molecule, effectively improving 



2 



the sampling by three order of magnitudes (there are 
1024 water molecules). Repeating the analysis with a 
longer equilibration time reproduced the same results. 
TIP5P data was collected up to 230 K, just before the 
approaching of the glass-transition [27?] . 



B. Density maximum 

The position of the maximum density was obtained 
from 1 ns long simulations after 10 ns of equilibration 
in the NPT ensemble. The temperature of maximum 
density was extracted by polynomial fitting around the 
maximum. Variations from the literature (see Table [I]) 
may be due to size effects and a different treatment of 
the electrostatics. The location of the TIP3P density 
maximum was obtained by running further simulations 
at lower temperatures. 



C. Energetics and tetrahedral order parameter 

The free energy of a configuration i is given by 



AFi = -k B Tlog(Pi), 



(1) 



where fee is the Boltzmann factor, T the temperature 
and Pi the population of the selected configuration. The 
enthalpy is estimated by summing up all pairwise con- 
tributions between the water molecules belonging to the 
same configuration. 

The tetrahedral order parameter for a water molecule 
i is calculated as 



j=i k=j+i 



COSIpjik 



-) 2 : 



(2) 



where j and k are any of the four nearest water molecules 
of i and ipjik is the angle formed by their oxygens [27]. 
The averaged value of this order parameter over an en- 
semble of water molecules is denoted as Q. 



III. RESULTS 
A. Structure forming capabilities 

Water structure forming capabilities were investigated 
by analyzing the hydrogen-bond network of each water 
molecule in the simulation box together with its first and 
second solvation shells. A maximum of four hydrogen- 
bonds per molecule was considered. A bond is formed 
when the distance between oxygens and the angle O-H- 
O is smaller than 3.5 A and 30 degrees, respectively [28] . 
Water structures were grouped into four archetypal con- 
figurations of population P^ : the fully coordinated first 



and second solvation shells for a total of 16 hydrogen- 
bonds (P4, see Fig. I] for a schematic representation); 
the fully coordinated first shell, in which one or more 
hydrogen bonds between the first and the second shells 
are missing or loops are formed (P4); the three coordi- 
nated water molecule (P3) and the rest (P2io)- Within 
this representation the sum over the four populations is 
equal to one for each temperature. 

In Fig. [2] the temperature dependence of the four mi- 
croscopic water structures is shown. Among the dif- 
ferent water models, the qualitative behavior is strik- 
ingly similar. Three main types of temperature scalings 
were observed: increasing population with decreasing 
temperature (enthalpically stabilized); increasing popu- 
lation with increasing temperature (entropically stabi- 
lized); with a maximum, where a turnover between en- 
thalpic and entropic stabilization takes place at a model 
dependent temperature. All four water configurations 
fall into one of these three main classes. The popula- 
tion of the fully ordered structure, P4, increases with 
decreasing temperature (Fig. [2j red empty circles) . Con- 
sequently, this configuration is enthalpically stabilized. 
This is not the case when defects in the hydrogen bond 
structure are introduced (P4, filled red circles). For this 
configuration the population increases with decreasing 
temperature until it reaches a maximum in correspon- 
dence to a rapid increase of the population of the fully- 
coordinated configuration. The maximum is located in a 
temperature range close to the temperature of maximum 
density of the model under consideration (dashed vertical 
line). Finally, both P3 and P210 are mainly entropically 
stabilized, showing larger populations at higher temper- 
atures. Taken together, these results indicate that spe- 
cific water configurations dominate at each temperature 
range: full-coordination extending to at least two solva- 
tion shells at low temperatures, four-coordinated config- 
urations with no spatial extension at intermediate tem- 
peratures and mainly disordered ones at higher temper- 
atures. 

Despite these similarities, an important difference 
among the models is the temperature range at which 
the relative configurations become dominant. For exam- 
ple, the maximum population of P| for the SPC model 
was observed around ~ 245 K (Table |TJ This is not the 
case for TIP4P/2005, where the maximum is located at 
a 40 K larger temperature. The same behavior was ob- 
served comparing the temperatures at which P4 and P3 
are equal (e.g., around 270 K for TIP4P/2005). These 
observations suggested that a temperature shift factor 
(AT S ) exists among the models. Using TIP4P/2005 as 
a reference, we found a temperature shift factor for each 
model ranging from 65 K to 6 K (see Fig. [3] and Ta- 
ble [I}. TIP4P/2005 was chosen as reference for its ability 
to reproduce the density curve {BJ. Applying this shift 
to the data allowed the superposition of all models onto 
four master curves, one for each structural configuration, 
as shown in the monochrome plot at the bottom right 
of Fig. [2] Our observation is consistent with previously 
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FIG. 2. Temperature dependence of water structure populations for seven classical water models. P4, P4, P3, and P210 are 
shown in red empty, filled red, blue, and light blue circles, respectively (see text for details). The gray stretch highlights the 
temperature difference between the calculated position of the temperature of maximum density (vertical dashed line, see also 
Table |T| and the experimental value at 277 K (solid line). The bottom right monochrome plot shows the superposition of all 
models after temperature shifting each data set (TIP4P/2005 data was used as reference). For each temperature, the sum over 
the 4 groups is equal to one. 
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FIG. 3. Structural temperature shift AT 3 with respect to the 
TIP4P/2005 model. TIP5P was excluded from the superpo- 
sition analysis (see text for details) . 



found phase diagram shifts among different water mod- 
els [221 [30] as well as in the presence of ions [3T] but in 
this case we could superimpose all models onto a master 
curve. Unfortunately, TIP5P had to be excluded from 
the superposition because all points show an increased 
curvature with respect to the other models, consistent 
with the increased curvature of the isobaric density at 1 
atm [8]. The structural temperature shift is larger for 
three-site models (Fig. [3]) with a spread of up to 65K for 
TIP3P. On the other hand, four-site models deviate less. 
Both SPC-E and TIP4P are characterized by a temper- 
ature shift with respect to TIP4P/2005 of around 20K. 
In general, models providing a better estimation of the 
position of the density maximum deviate less. 

In order to check the robustness of the Pi overlap with 
the hydrogen bond definition, the recent definition of 
Skinner and collaborators [32] was applied. Fig. [4] shows 
that the overlap between the curves is independent from 
the hydrogen bond definition. Moreover, the tempera- 
ture shifts calculated in this case are very similar to the 
ones reported in Table [T] It is interesting to note that 
this bond definition is very different from the one used in 
Fig. [2] being based on an empirical fit of the electronic 
occupancy 32 . 



B. Stabilization of the fully coordinated 
configuration 

At all temperatures, water models with smaller shifts 
provide an improved stabilization of the fully coordi- 
nated configuration. (Alternatively, it can be said that 
these models destabilize poorly hydrogen-bonded config- 
urations). To make this point clearer, the free energy 



of the fully coordinated configuration at 230K was cal- 
culated (Fig. [5^, and Methods). At this temperature P4 
is appreciably large for all water models. Comparison 
with the temperature shifts of Fig. [3] indicates a remark- 
able correlation where even the small differences between 
SPC-E and TIP4P are respected (Pearson correlation co- 
efficient of 0.99). The correlation decreases when consid- 
ering only the enthalpy, as shown in Fig. [SJd (correlation 
of 0.94, see Methods). Interestingly, enthalpy and free 
energy correlate very well within the same model family. 
This is particularly clear when looking at the three sites 
models (i.e., the trend for TIP3P, SPC and SPC/E), sug- 
gesting a different cntropic contribution between three 
and four sites models which is systematic. 

Finally, the average value of the tetrahedral order pa- 
rameter Q [27] of the fully coordinated configuration cal- 
culated at the same temperature is shown in Fig. (5ji. In 
first approximation, the parameter correlates well with 
the structural shift although not as good as the free en- 
ergy (correlation coefficients of 0.98). 



C. Comparison with the position of the density 
maximum 

It is worth commenting on the relation between the 
structural temperature shifts found in this work and the 
mo del- dependent temperature of maximum density. As 
shown in Fig. [6] and Table [I] the relationship between the 
structural AT S and the density ATrf eras jj y temperature 
shifts is linear within the three or the four-sites models 
(filled and empty circles in Fig. [6|. However, when com- 
paring all models together using TIP4P/2005 as reference 
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FIG. 4. Overlap of the Pi data when a Skinner definition [32] 
for the hydrogen bond is used. 
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FIG. 5. Comparison of water models with respect to the 
fully coordinated configuration at 230 K. (a) The value of the 
free energy, (b) Average enthalpy, (c) Average value of the 
tetrahedral order parameter Q. The value of Q for the case 
P4 = 0.2 is shown as gray filled circles. 



a small systematic deviation is observed (filled circles and 
crosses in the figure) . This is due to the relation that ex- 
ists between the populations Pi and the density. It is 
noted that the relative position of the P4 maximum with 
respect to the temperature of maximum density (dashed 
line in Fig. [2j see also Table [I]) depends on the model 
family. For four-sites models the two temperatures are 
identical, while for three and five-sites models the maxi- 
mum of P| is found at a higher and a lower temperature, 
respectively. This behavior might be connected with the 
systematic deviations between free energy and enthalpy 
for the different water models (Fig. [5^,-b) . 

To better elucidate the connection with the tempera- 
ture shifts, the radial distribution function (RDF) was 
calculated. At 270 K the RDF for the various models 
shows a different structural signature (Fig. [7^,). As ex- 
pected, only the curves for TIP4P/2005 and TIP4P-Ew 
overlap, being the two models very similar. In Fig.[7j3 the 
RDF was recalculated at temperatures shifted according 
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FIG. 6. Comparison between the structural temperature 
shifts (AT 3 ) and the position of the density maximum 
(ATdensity) ■ Four-site models were compared to TIP4P/2005 
(filled circles). Three-site models were compared to SPC/E 
(empty circles). Crosses refer to the case when TIP4P/2005 
was used as reference for the three-site models. 



to AT S using as reference TIP4P/2005 at 270 K. The fig- 
ure shows that all models with the same geometry per- 
fectly overlap (e.g. TIP4P, TIP4P-Ew and TIP4P/2005, 
blue lines) , suggesting that model reparametrizations act 
as an effective shift in temperature space. On the other 
hand, changes in the geometry or the number of sites af- 
fect the general shape of the radial distribution function 
and consequently the density. 

Similar conclusions can be deduced when calculating 
the tetrahedral order parameter Q for temperatures at 
which the fully coordinated structure has the same prob- 
ability for all models (P4 — 0.2, gray filled circles in 
Fig. [5b). Q takes the same value within a given fam- 
ily but it is influenced by the change of the molecular 
geometry, indicating that the structure corresponding to 
fully coordinated waters depends on the model family. 



IV. CONCLUSIONS 

In conclusion, we found that seven among the most 
used classical water models are characterized by very sim- 
ilar hydrogen-bond structure-forming capabilities up to a 
temperature shift. All models but TIP5P perfectly over- 
lap onto a master curve when this shift is applied. This 
behavior does not depend on the hydrogen-bond defini- 
tion. Our findings suggest that model reparametrization 
acts as an effective shift in temperature space. On the 
other hand, changes in the geometry or the number of 
sites cannot be fully reconducted to temperature shifts 
alone as shown by the analysis of the density as well 
as the radial distribution function. As such, although 
the hydrogen bond topology is universal when applying 
a certain temperature shift, this is not the case for the 
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FIG. 7. Radial distribution function (RDF) (a) at 270 K; (b) after the application of the temperature shift AT S (see Table I), 
taking the TIP4P/2005 model at 270 K as reference. TIP3P, SPC and TIP4P families are shown in red, light blue and blue, 
respectively. 



structure, each model family being characterized by its 
own signature. 

We found that the three water models optimized to 
reproduce the position of the density maximum (i.e., 
TIP4P-Ew, TIP4P/2005 and TIP5P) systematically im- 
prove the stabilization of fully coordinated water con- 
figurations with an extension of at least two solvation 
shells. Based on this observation, we speculate that the 
improvements of these models for biomolecular simula- 
tions [131 [TS] are connected to the higher stabilization of 
ordered water. This property has important implications 
for the solvation of biomolecules, changing the balance 



between solute-solute and solute-solvent interactions. 

Development of improved force-fields strongly depends 
on this balance. Our analysis provides a microscopic and 
reductionist approach to face this challenge. 
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TABLE I. Temperature of maximum density calculated from our simulations (TMD), as found in the literature (TMD re '), the 
structural temperature shift (AT S ) and the temperature at which P4 is maximum for the seven water models investigated in 



this work. 

Water model TMD TMD re/ ATs T max( p 4 ') 

TIP3P 199 1820 65 229 

SPC 226 2280 42 247 

SPC/E 250 241 [33] 18 275 

TIP4P 256 248[II] 20 268 

TIP4P/2005 280 2780 287 

TIP4P-Ew 273 274[l8] 6 281 

TIP5P 282 285[3J] n.a, 269 



